# -*- coding: utf-8 -*-
from __future__ import unicode_literals

"""
Created on Tue Dec 13 21:28:59 2016

@author: lihan
企业状态分析 海关数据
"""
import numpy as np
import pandas as pd
from scipy.optimize import leastsq
import pylab as pl


def fitfunc(p, t):
    return p[0] / (1 + p[2] * np.exp(-p[1] * t))


def lossfunc(p, y, t):
    return fitfunc(p, t) - y


def lifecurve(t_data, dd):
    aa = np.array(dd)
    data = pd.DataFrame(aa)
    # data = pd.read_excel(io=filename, header=None)  # 读入海关数据，
    Y = np.array(data, dtype=float)  # Y 是海关scale，需按月份排列，需剔除不稳定的跳跃点
    # 取前三个序列最小二乘求logistic曲线的参数[N，alpha，r]的初始值
    t = np.array([1, 2, 3])  # t为月份时间序列
    Y = Y[Y > 10]
    YY = Y[0:3]
    if len(Y) < 3:
        return 0, 0
    N = float((2 * YY[0] * YY[1] * YY[2] - (YY[1] ** 2) * (YY[0] + YY[2])) / (YY[0] * YY[2] - YY[1] ** 2))  # scale饱和值
    y = np.log(YY) - np.log(abs(N - YY))
    a = abs(np.polyfit(-t, y, 1))
    r = a[0]
    alpha = np.exp(a[1])
    # 非线性回归求求logistic曲线的参数[N，alpha，r]的最优值
    t = range(len(Y))
    t = np.array(t)
    Y = np.array(Y)
    beta, res = leastsq(lossfunc, np.array([N, r, alpha]), args=(Y, t))
    N = beta[0]
    r = abs(beta[1])
    alpha = abs(beta[2])

    T0 = np.log(alpha) / r  # T0是成长期前期与成长期后期的分界点
    T1 = (np.log(alpha) - np.log(2 + np.sqrt(3))) / r  # T1是初创期与成长期前期的分界点
    T2 = (np.log(alpha) + np.log(2 + np.sqrt(3))) / r  # T2是成长期后期与成熟稳定期的分界点
    if np.isnan(T0) or np.isnan(T1) or np.isnan(T2):
        return 0, 0
    if T1 > t_data:
        pass
    elif T1 <= t_data:
        k = T1 / t_data + 1
        T0 = T0 / k
        T1 = T1 / k
        T2 = T2 / k
    else:
        return 0, 0
    x = pl.frange(t_data, 240 + t_data, 12)  # x时间序列
    yy = N / (1 + alpha * np.exp(-r * x))  # 回归拟合的scale
    yyy = yy * r * (1 - yy / N)  # scale的增长率

    fenge = [T1, T0, T2]  # 分界点时间，对应的实际月份是T+基期（样本最小的月份如14年1月）
    x = map(lambda z:z/12, x)
    fenge = map(lambda z:z/12, fenge)
    coordinate = [x, yyy]
    return coordinate, fenge

# if __name__ == '__main__':
#     aa = [
#         627062400.0,
#         658598400.0,
#         690134400.0,
#         721756800.0,
#         753292800.0,
#         784828800.0,
#         816364800.0,
#         847987200.0,
#         879523200.0,
#         911059200.0,
#         942595200.0,
#         974217600.0,
#         1005753600.0,
#         1037289600.0,
#         1068825600.0,
#         1100448000.0,
#         1131984000.0,
#         1163520000.0,
#         1195056000.0,
#         1226678400.0,
#         1258214400.0,
#         1289750400.0,
#         1321286400.0,
#         1352908800.0,
#         1384444800.0,
#         1415980800.0]
#     print lifecurve(-48, aa)
